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Introduction 


Upwind differencing, while trivial for a diagonalized hyperbolic system, 
is difficult to achieve when the difference scheme has to be written in 
conservation form. The oldest and most complicated version is due to 
Godunov [1] [2]; the increasing popularity of upwind differencing recently 
has lead to a variety of simpler implementation techniques. A review of 
these is given by Harten, Lax and Van Leer [3], 

Among the recent additions to the family of upwind conservative schemes 
the method of Engquist and Osher [4] [5] is closest to the original Godunov 
scheme, while the method of Roe [6] [7] offers the greatest simplification 
In the present paper the differences between these schemes are discussed on 
the basis of the inviscid Burgers equation. The first-order accurate schemes 
are explained in Sections 2, 3, and 4 for the homogeneous equation; Section 6 
describes how to include a source term and Section 7 how to achieve second- 
order accuracy in a two-step format. Their rendition of a stationary shock 
and a transonic expansion is discussed in Section 5 and illustrated in the 
numerical experiments of Section 8. Section 9 rounds off with recommendations 
regarding the application of these schemes to single conservation laws and 
systems of conservation laws. 

2. Godunovas method 

Godunov’s [1] [2] method for integrating a hyperbolic system of conser- 
vation laws 

(1) + [f(u)]^ = 0, 

is a scheme in conservation form: 


(u^"^^-u^)/At + {F(u^,u^^^) - F(u^^^,u^)}/Ax = 0; 


( 2 ) 
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here represents the average value at time t^ = nAt in the computational 

zone centered on = iAx. The numerical flux-function 

Godunov scheme is taken to be the flux value arising at x^^ ^ i in the exact 
solution of the initial-value problem with a piecewise uniform initial distri- 
bution 

(3.1) u^(x) = u^ x^ - Ax/2 < X < x^ + Ax/2. 

That is, if 


(3.2) 

u(x,t) 

= v(x/t;u^,Uj^) , 

is the 

(weak) similarity solution 

of the Riemann problem with initial 

values 




("L 

X < 0 

(3.3) 

u = 



K 

X > 0, 


then 


(4) 




For Burgers* equation in the inviscid limit, 

(5) Uj. + = 0, 

the similarity solution to the initial-value problem (3.2) is 


( 6 . 1 ) 


(expansion) 


V = u^ 

V = x/t 

V = Uj^ 


x/t < u^ 

n < x/t < u 
L 

x/t > u^ 


( 6 . 2 ) 


(shock) 




V = u„ 


x/t < Ug = j(Uj^ + Uj^) 
x/t > Ug 
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For numerical reasons it is useful to distinguish in the formula 

for 

(i) fully supersonic/subsonic 


ll) three cases: 


(7.1) 




■luj 


0 



0 , 


(ii) transonic expansion 


(7.2) 




< 0 < Ur , 


(iii) transonic shock 



f 4 ut 

u, > u - > 0 > u, 

(7.3) 


L S — — 1 


iiul 

\> 0 >^s ^ 


These cases are illustrated in Figures 1, 2 and 3. The references to a 
sonic speed ( = 0) arises from the use of (5) in transonic aerodynamics. 
In order to combine the formulas (7) in one compact algorithm, we 


follow Engquist 

and Osher by introducing 

(8.1) 

u^ = max(u,0) , 

(8.2) 

u = min(u,0) , 

(8.3) 

+ 

u = u + u , 

(8.4) 


|u| = u - u ; 

we then have 



Fg(Uj^.Ur) = max [4(u^)^ , 4(Uj^)^] . 


(9) 
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3. The Engquist-Osher scheme 

The Engquist-Osher [4] [5] scheme for integrating (1) also has the 
conservation form (2) and tacitly assumes the initial-value distribution 
to be (3.1); its numerical flux-function is 


(10) - f(u^)- / A (u)du 

U 

= i[f (Uj^) + f (Uj^)] -if ^|A(u)|du, 

where the matrices A^(u) , A~(u) and |a(u)| are related to 

(11) A(u)= df/du 


by an extension of (8) . For precise definitions and a description of the 
integration path used for the integrals in (10), see [5] or the 
review [3]. 

For Burgers* equation (5) the above recipe boils down to 


( 12 ) 


= 4u^ + / ^ u du 


1 2 +, 

= iu - j u du 


1 

2 


(4uJ+4u^) - 



The integrals over u in (12) are defined in phase space, without 
reference to any mapping onto the (x,t) plane. We may, however, intro- 
duce a mapping in the style of (3.2), namely. 
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(13) u(x,t) = w(x/t;u^,Uj^) min(u^,u^) £ u 5 max(Uj^,Uj^) . 

and consider w(x/t;u ,u ) an approximation to the exact solution 

L K 

v(x/t;u^,u ) of the Riemann problem (3.3). It follows that w is the 
L K 

following function of x/t : 


(14) 


W = x/t < , 

w = x/t min(Uj^,Uj^) < x/t < max(Uj^,Uj^) , 

w = Ur x/t > Ur . 


In relating Frq(Uj^.Ur) to w(0;Uj^,Ur) after the example of (4), 

caution is required, since, for u >u , w becomes multi-valued in the 

L K 

domain u^ > x/t > u„ . Specifically, we have three branches 
L — K 

(15) w^^^*u^, w^^^=x/t, w^^^=Uj^ ^ x/t £ 


The picture associated with this case is that of an overturned centered 
compression wave or folded characteristic field (see Figures 4 and 5); 
in the exact Riemann solution ( 6 ) such a wave would be replaced by a shock 
discontinuity. 

The proper formula for ^eq^^L’^R^ terms of w(x/t;Uj^,Uj^) , 

equivalent to ( 12 ), io 

(16) fEO<"L'V 

k 

where the sum is taken over all branches present at x/ 1 = 0 . 
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The difference between the Godunov and Engquist-Osher schemes lies 

entirely in the treatment of a transonic compression (iii) . The latter 

scheme is the simpler one, since it does away with one test, namely, a 

+ 2 - 2 

test for the sign of or for the maximum of and 


4. Roe^s method 

Roe’s [6] [7] method for integrating (1) again has the conservation 
form (2) and employs the initial-value distribution (3.1), Its numerical 
flux-function is defined as 

(19) " f [w(0;u",u^^^)l , 

where w(x/t;u^,u^) is the exact solution of the Riemann problem (3.3) 
for the locally linearized system 

(20.1) w^ + A(u^,Uj^)w^ = 0. 

The approximate Jacobian A(u ,u ) is constructed such as to satisfy 

L R 

the discrete version of (11) 
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( 20 . 2 ) ’ 

For a scalar equation like (5), eq. (20.2) uniquely determines 
A(Ul,u^): 

(21) A(uj^,Uj^) = = Ug; 

thus, the Riemann problem for the linear equation (20.1) with (21) has 
the similarity solution 



The difference with Godunov’s method lies in the treatment of 
an expansion: where the exact Riemann solution, used in Godunov’s method, 
would include an expansion fan, Roe’s method puts in a so-called expan- 
sion shock (see Figure 6) • 

The numerical flux-function (23) deviates from the Godunov flux (7) 
only in case (ii) , when the expansion is transonic. The computational 
simplification is, again, the elimination of one test. 

Rewriting (23) as an approximation to the Enquist-Osher flux (12): 


we see that, for Burgers’ equation. Roe’s scheme is identical to the scheme 
of Murman and Cole [8], used in transonic aerodynamics. The latter is 
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known to occasionally yield numerical results that include a (physically 
inadmissible) expansion shock. This is a direct consequence of the ad- 
mission of an expansion shock in the underlying approximate Riemann 
solution. 


5. Steady-shock and sonic-point representation 

The fluxes in the schemes of Godunov and Engquist-Osher differ only 
on meshes where the data constitute a transonic compression; therefore, 
numerical results from these schemes differ only if a transonic shock, 
in particular, a stationary shock, is present. Likewise, the results of 
Roe’s scheme will differ from those of Godunov’s scheme only if a transonic 
expansion is present. 

For Burgers’ equation (5), Godunov’s scheme admits the following 
stationary discrete representation of a stationary shock connecting the 
states u^ > 0 and u^^ = -u^^ <0; 


= “l 


i < -1, 


(25) 


“O = ^ 




u. = -Ut i > 1. 

This is the only stationary distribution admitted by the scheme. 

The interior value u^^ indicates the sub-grid position Xg of the 
shock; by conservation we have, (see Figure 7a): 


(26.'1) iY(Xg + 4Ax) + u^(iAx-Xg) = tij^Ax, 

or 


(26.2) 


Xg = iAx 



- 12 - 



Flgure 7. Fitting a shock discontinuity to stationary discrete shocks 
yielded by the schemes of Godunov or Roe (a) and Engquist- 
Osher (b) . The numerical values u^ and u^ represent the 
average value of u in the intervals (-JAx,jAx) and (-jAx, 
-|Ax) . A sub-grid distribution representing a shock transition 
from Uj^ to at Xg must have the same integral as the 

numerical solution. In case (a) the integrals from -jAx to 
jAx are compared, in case (b) from -JAx to -|Ax. 
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A shock standing exactly on the boundary of a zone will be represented 
without any interior value. 

The Engquist-Osher scheme admits discrete steady shock-profiles 


with. 

in general. 

two Interior 

states ; 




u . 
1 


1 < -1, 







0, 

(27.1) 










0 > > • 






i > 2. 


where 


are constrained by 


(27.2) 







Again, this is the only stationary distribution admitted by the scheme. 

The pair of interior states indicates the sub-grid position of the 
shock; by conservation we have (see Figure 7b) : 

(28.1) Uj^(Xg + 4Ax) + Ax-Xg) = (Uj^+Uj^)Ax, 

or 

(28.2) Xg = 4Ax(Uj^+Uj^+Uj^)/uj^. 

A shock standing exactly in the middle of a zone will be represented 
with only one interior value. 

The representation of a stationary shock by Roe's scheme is the same 
as by Godunov's scheme, at least for Burgers' equation. However, Roe's 
scheme also admits as a stationary solution the expansion shock 
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u. = u < 0 i < -1, 

X Li 

(29) 

= -Uj^ >0 i 

The tendency of the scheme to produce such a discontinuity spoils the 
smoothness of the numerical solution near a sonic point. The other 
schemes do not have this problem. It is easily solved, anyway, by intro- 
ducing in the approximate Riemann solution (22) some spreading of the 
expansion shock- For Burgers’ equation, Roe's scheme essentially will turn 
into Godunov’s scheme; for a nonlinear hyperbolic system, however, it will 
still remain a simpler scheme. 

A transonic shock profile, steady or not, obtained with any of the 
upwind schemes, has the property that the interior zones can not influ- 
ence the exterior solution. Inversely, in a transonic expansion computed 
with the Roe-Murman-Cole scheme, the zone with the value closest to the 
sonic value can not be Influenced by the rest of the grid (see Figure 8). 

6- Inclusion of a source term 

When approximating the equation 

(30) u^ + [f(u)]^ = s(x) 

with any upwind scheme, it is not sufficient to add the source term to 
the right-hand side of (2) ; we must also change the initial-value distri- 
bution implied in our numerical model- Instead of assuming that it is 
piecewise uniform, as in (3.1), we rather take it to be piecewise 
stationary , that is, 

(31.1) [f(u^(x))]^ = s(x) x^ - Ax/2 < X < x^ + Ax/2, 
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/VI R 


Figure 8, Two (x,t) diagrams showing the influence, through Godunov’s 
scheme, of the adjacent zones on a zone inside a shock (a) 
and the influence, through Roe’s scheme, of an almost sonic 
zone on the adjacent zones (b) . In (a), the fluxes at the 
boundaries of zone M do not depend on u^^; in (b) they 
depend only on u^. 
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U 



Figure 9. Stationary solutions of eq. (30) for some source term 

s(x). In each zone the heavy line traces the positive and 
the negative branch of the smallest continuous steady solu- 
tion. The average value on these branches is indicated by 
the upper and lower boundary of the shaded area. If the 
average value in a zone falls in the range of the shaded 
area, a continuous stationary solution can not be realized. 
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with 

(31.2) 



+ Ax/2 


n 


u 


(x)dx 


x^ - Ax/2 



This extension preserves a fundamental property of the homogeneous scheme, 
namely, that a zone-average can change only through the finite-amplitude 
waves entering from the zone boundaries. Moreover, if the numerical solu- 
tion globally tends to a steady state, it will be the zone-averaged exact 
steady state almost everywhere (that is, provided that the scheme, like 
the three considered here, can render a shock transition in a finite 
number of zones) - 

As before, the initial values on either side of a zone boundary become 
the arguments of the numerical flux-function. The full scheme reads 



under the constraint (31.2). This constraint, however, is not strong enough 
to define a unique distribution in each zone, since the distribution may 
(and, in some zones, must) contain a discontinuity. Uniqueness can be 
achieved by .selecting the distribution with, say, the weakest possible shock. 



Figure 10. Examples of steady structure with a discontinuity, in 

the zones shown in Figure 9. A (steady) shock connects 
the upper and lower branch of the smallest continuous 
steady solution for the zone considered. In zone j 
two different possibilities are shown. 
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If we take Ax small enough to ensure that s(x) does not change 
sign more than once in any zone, the following algorithm for the piecewise 
construction of u^(x) is adequate. First, determine for the zone under 
consideration the smallest continuous solution (see Figure 9). Let 
be the average value on the positive branch of this solution; then the 
stationary distribution inside zone i will be continuous if 


> 


If this condition is satisfied, we search, iteratively, for the continuous 
distribution with the proper zone average. If it is violated, the sta- 
tionary distribution will include parts of the upper and lower branch 
of the smallest continuous distribution, connected by a shock positioned 
so as to achieve the proper average value (see Figure 10) . 

When using the upwind scheme to approach a globally stationary solution 
of (33), any zone containing a sonic point must be treated with special care, 
since the chance of numerically realizing the exact transonic structure 
without a shock is zero. Specifically, in order to prevent the zone- 
boundary values from flipping sign, making global convergence impossible, 
smoothing may be introduced (see Figure 11) . 

In a zone containing a stationary shock, no particular interior structure 
is needed to ensure global convergence, since the zone cannot influence its 
neighbors (see earlier Figure 8a) . After convergence one may insert the 
proper structure by enforcing continuity across the zone boundaries (see 
Figure 12). 

In the numerical experiment of Section 8 we approximated the above 
procedure, in allowing only s^ to enter the zone-boundary values. 
Specifically, we used 



Figure 11. (a) A sequence of stationary structures close to the transonic expansion (graph IV), ordered 

by decreasing zone-average. Going from III to IV, and from IV to V, the slight change in 
zone-average causes one boundary value to flip its sign. (b) The solutions have been smoothed 
on the side where the shock occurs, so that the boundary values now vary continuously with the 
zone-average. The structures are not stationary but will allow convergence to the smooth 
transonic solution IV. 
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Figure 12. Shock structure in a globally stationary solution. Global 
convergence of the zone averages to a steady solution con- 
taining a shock has been obtained, say, with Godunov* s 
scheme. For zone k containing the shock the scheme had 
adopted the structure (a), based on the smallest continuous 
steady solution for that zone. This may now be replaced 
by the stationary solution (b) that smoothly connects to 
the solution in zones k-1 and k+1. 
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(36) sgn(u”) / (u”)^±min{|s^|Ax, (u^)^}sgn(s^) • 

In a zone that includes a sonic point the error in the boundary values 
is 0(Ax); elsewhere in the smooth part of the numerical solution it is 
0((Ax)2). 

The technique of computing the boundary values in a zone from a 
locally stationary solution can easily be extended to systems of conserva- 
tion laws of the form 

(37) u^ + [f(u,x)]^ = s(x). 


7. Second-order upwind schemes 

Any first-order upwind scheme for Eq. (30) can be changed into a second- 
order method by first advancing the cell-boundary values, to be used in the 
numerical flux-function, to the intermediate time level t^"^ = + iAt. 

In obtaining these values, the interaction between cells can be fully ignored. 
This observation, due to Hancock [9], has led to a drastic simplification of 
second-order upwind schemes since these first were formulated for systems of 
conservation laws by Van Leer [10]. 

We assume that the initial values form a piecewise linear distribution : 


(38) 


(5u)” 

u (x) = u. + (x-x.) 


X i < X < X ^ , 


with 


j, r. ^n . n n n n . 

(6u)^ = ave(u^^^-u^, u^-u^__^); 


( 39 ) 
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ave(a,b) is an averaging procedure to be specified later. We particularly 
need the initial boundary values inside cell i 


(40.1) 




These boundary values are now advanced to t 


n+^ 


by 


(40.2) 


n+4 


At 


At 


(i±4)+ '^(i±i)+ 2Ax 2 


The full scheme becomes 


(41) (u"‘^^-u")/At 


( / n-fi n+J \ / n*fj n+i \ ) 

j^(''(l+i)-’''(i+i)+) ^(“(i-i)-*'‘(i-4)+)j 


/Ax = s 


i“ 


The function ave(a,b) is chosen such that it tends toi 4(^+b) if 
a and b are subsequent finite differences of a smooth solution, but tends 
to the smallest value where the solution is not smooth. Examples can be 
found in [11], [10]; we chose a refined formula due to Van Albada [12] 


(42.1) 


ave(a,b) = 


2 2 2 2 
(b + c )a + (a + c )b 

2 2 2 
a"^ + b + 2c 


2 3 

where c is a small bias of the order 0((Ax) ). 

The weighted averaging prevents central differencing across a discon- 
tinuity in the solution or in its first derivative which would lead to 
numerical oscillations. It is an effective way of administering artificial 
dissipation wherever needed and nowhere else. By rewriting (42.1) as 


(42.2) 


ave 


- a+ b (a - b) 

(a,b) — n “ o o 


2 2 2 
a + b + 2c 


we see that, wherever the solution is smooth, the artificial-viscosity 

2 

coefficient generally is of the order 0((Ax) ). In a smooth extremum the 

2 

coefficient grows to 0(Ax); the bias 2c in the denominator prevents a 
further increase of the viscosity that could lead to an undesirable clipping 
phenomenon (see [13, Appendix]), 
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With regard to Burgers’ equation, an acceptable expression for the 
bias is 


(43) 


c^ = (u -u , )^(Ax)^/(x -X . )^, 
max raxn' ^ ^ min 


where u and u are certain upper and lower bounds of the numerical 
max min ^ ^ 

solution, fixed a priori or determi 
is the length-scale of the problem. 


solution, fixed a priori or determined at each time level, and x -x . 

max min 


8 , Numerical comparison 

The performance of the three schemes was tested on the basis of the 
periodic initial-value problem 

(44.1) + (4u^) = (Tr/2)sin[27T(x- C)] 0<x<l, 0<5 « 1 

t X 

(44.2) u(x,0) = 0, 

(44.3) u(0,t) = u(l,t). 

The solution tends to the steady state 

j+sin7T(x-5) 0^x< 5 + i 

(45) u(x,“) = I 

(-sin7T(x-5) 5 + 

including a sonic point at x = 5 and a shock at x = x^ = ^ + For 
Ax we chose a value of 1/16; the zones were centered on x^ = (i- J)Ax, 
i - 1,...,16. For C we chose 0., ^Ax, and jAx, in order to achieve 
different sub-grid shock positions. The Courant-Friedrichs-Lewy condi- 
tion on the timestep, based on the maximum characterisitc speed (=1) in 


the steady solution is 
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(46) ^ < i; 

accordingly we used At = -jAx or At Ax. 

Table 1 shows the number of steps it took the schemes to con- 

verge according to the L^- criterion 

16 

('^•7) I |u“-u^"^| < 0, 

i=l 

—3 —6 

with e = 1 X 10 or 1 x lo , and At = jAx or (for Godunov *s scheme 
only) At = Ax. When the time-step is doubled, appears to be halved. 

The Lj^-error was evaluated with respect to the zone-averaged exact 

solution, for the smallest value of e. The sub-grid shock positions in 
Table 1 were calculated with aid of eqs. (26.1) and (28.1), approximately 
valid because the source term is small in the shock region. 

The converged solutions for e = 1 x 10 ^ are listed in Table 2; 
these are independent of At. The zone-averaged exact solution is given 
for comparison. For ^ = JAx and JAx the error in the zone with the sonic 
point is comparable to the value in the zone, as anticipated with the use of 
Eq. (36). 

The results of Roe’s scheme happen to be identical to those of Godunov’s 
scheme; in particular, no expansion shock shows up. The reason is that, 
with the use of the boundary values (36) , a sonic point in a zone is always 
moved to a boundary; in this case Roe’s flux function (23) yields the same 
value (=0) as Godunov’s (7.2). Furthermore, the initial average value in the 

zone ultimately containing the sonic point was close to the steady value 
to begin with. 

The results of the Engquist-Osher scheme are identical to the Godunov- 
Roe results except for the expected spread in the shock profile for 5=0 
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E 


n 

5 /Ax 

e = 1 X 

10 ^ 

E = 1 X 10 ^ 

e = 1 : 

s 10"^ 



G, R 

EO 

G, R 

EO 

G, R 

EO 

G, R 

EO 

0 

62 

61 

112 [55] 

111 

8.8 X 10"^ 

4.6 X 10 ^ 

0 

0 

i 

68 

66 

138 [70] 

136 

9.6 X 10"^ 

1.8 X 10“^ 

.23 

.26 


52 

52 

88 [42] 

88 

4.6 X 10 ^ 

4.6 X 10 ^ 

■J 

1 

2 


Table !• Number of steps till convergence, L^^-error E^, and 

steady-shock position Xg , for the initial-value problem 
(44), solved with the first order schemes of Godunov (G) , 
Engquist-Osher (EO) and Roe (R) . Mesh; Ax = 1/16, 

At = jAx or Ax (G only, numbers between brackets). 



















.13795 

.30373 

.47702 


.13795 

.30373 

.47702 


09706 

28902 

47064 


.09778 

.25516 

.43185 


.09778 

.25516 

.43184 


.63587 

.77180 

.87889 

.95276 

.99044 

-.99044 

-.95276 

-.87889 

-.77180 

-.63587 

-.47702 

-.30373 

-.13795 


.63587 

.77180 

.87889 

.95276 


rarjCiW 


-.69352 

-.95276 

-.87889 

-.77180 

-.63587 

-.47702 

-.30373 

-.13795 



.59602 

.73869 

.85367 

.93631 

.88033 

-.42699 

-. 96441 

-.89927 

-.79997 

-.67047 

-.51616 

-.34426 

-.16827 



5 /Ax = 4 


exact 

G, R, EO 

exact 

.04899 

0.00000 

0.00000 

.24259 

.19321 

.19478 

.42687 

.37899 

.38207 

.59474 

.55021 

.55468 

.73976 

.70028 

.70597 

.85635 

.82344 

.83013 

.94003 

.91496 

.92240 

.98759 

.97132 

.97921 

-.49739 

.00000 

-.00000 

-.96847 

-.97132 

-.97921 

-.90254 

-.91496 

-.92240 

-.00192 

-.82344 

-.83013 

-.67048 

-.70028 

-.70597 

-.51328 

-.55021 

-.55468 

-.33635 

-.37899 

-.38207 

-.14649 

-.19321 

-.19478 


±c solutions for the experiments 
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and jAx. The sub-grid shock positions are slightly, but no significantly, 
more accurate than for the other schemes. Likewise, convergence with 
Engquist-Osher scheme is faster, but not significantly so. 

Better results (Tables 3 and 4) were obtained with the second-order 
versions of the schemes. Convergence was achieved in 10-35% fewer steps 
than with the first-order schemes, and the local error near the sonic point 
is now of the same order as elsewhere outside the shock. 

The Engquist-Osher scheme still yields a shock structure with two 
interior zones. For the second-order Godunov scheme we checked that the 
use of algebraic averaging of differences 

(48) ave(a,b) = 4(a+b), 


causes the numerical shock to overshooot and undershoot. 

The results of Roe*s scheme are practically, but not exactly, identi- 
cal to the Godunov results, indicating a slightly different transient 
behavior. Again, the internal structure (36) in the zones, in combination 
with the particular choice of initial values, provides a safeguard against 
the occurrence of an expansion shock. 

In order to make the problem more challenging, we changed the initial 
values (44.2), for ^ = 0 ^ into 


(49) 


Li 


i = 1, ... ,8 
i = 9, — ,16, 


including an expansion shock at x = 0. The results are listed in Tables 5 
and 6. Among the first-order schemes Roe’s scheme now requires 30% more 
steps than the other schemes: apparently, the expansion shock is not so 

easily dissipated by Roe’s scheme. Among the second-order schemes the dis- 
crepancy gets worse: the results of Roe’s scheme quickly converge to the 

wrong solution, including the full initial expansion shock. 
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N 

e 





- 

i)/Ax 



G2,R2 

G2a 

E02 

G2, R2 

G2a 

E02 

G2, R2 

G2a 

E02 

0 

75 

76 

75 

1.4 X 10"^ 

3.1 X 10"^ 

2.9 X 10"^ 


0 

0 

4 

89 

92 

89 

1.3 X 10~^ 

2.2 X 10“^ 

4.7 X lO"^ 

.24 

.19 

.26 

i 

79 

79 

79 

1.3 X 10"^ 

1.7 X 10"^ 

1.3 X 10“^ 

i 

i 

i 


Table 3. The quantities N^, and Xg (e = 1 x 10~^) for the initial-value 
problem (44), solved with the second-order two-step schemes based 
on the numerical flux-functions of Godunov (G2, G2a) , Engquist- 
Osher (E02) and Roe (R2), In G2a the algebraic average (48) is 
used instead of (42.1). Mesh: Ax = 1/16, t = iAx. 
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00 

i 

> 

II 
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G2 

G2a 

R 

EO 

exact 

1 

.09814 

.09814 

.09815 

.09814 

.09786 

2 

.28976 

.28967 

.28976 

. 28975 

.28982 

3 


.46981 

.47012 

.47012 

.47064 

4 

.63242 

.63201 

.63242 

.63239 

.63337 

5 

.77076 

.76909 

.77076 

.77086 

.77177 

& 

.87879 

.88201 

.87880 

.87839 

.88851 

7 

.95647 

.92357 

.95648 

.95804 

.95540 

8 

.98761 

1.20652 

.98761 

.76751 

.99359 

9 

-.98761 

-1.20652 

-.98761 

-.76751 

-.99359 

10 

-.95647 

-.92357 

-.95648 

-.95804 

-.95540 

11 

-.87879 

-.88201 

-.87880 

-.87839 

-.88851 

12 

-.77076 

-.76909 

-.77076 

-.77086 

-.77177 

13 

-.63242 

-.63201 

-.63242 

-.63239 

-.63337 

14 

-.47012 

-.46981 

-.47012 

-.47012 

-.47064 

15 

-.28976 

-.28967 

-.28976 

-.28975 

-.28982 

16 

-.09814 

-.09814 

-.09815 

-.09814 

-.09786 


Table 4. 


Converged numerical solutions and zone-averaged asymptotic 
solutions for the experiments of Table 3, with G = 1 x 10 
At = -jAx. 


















N 

e 

Ui , 

OO 


i 



5/Ax = 

i 


5/Ax = i 


G2 

G2a 

R 

EO 

Exact 

G2, R2, EO 

G2a 

exact 

1 

.04916 

.04917 

.04916 

.04916 

.04899 

0.00000 

0.00000 

0.00000 

2 

.24249 

.24244 

.24249 

.24248 

.24259 

.19468 

.19465 

.19478 

3 

.42635 

.42610 

. 42635 

.42639 

.42687 

.38160 

.38141 

.38207 

4 

.59379 

.59344 

.59379 

.59364 

.59474 

.55379 

.55348 

.55468 

5 

.73871 

.73726 

.73871 

.73933 

.73976 

.70492 

.7^0373 

.70597 

6 

.85462 

.85711 

.85462 

.85231 

.85635 

.82844 

.83000 

.83013 

7 

.94080 

.91445 

.94080 

.95106 

.94003 

.92282 

.90418 

.92240 

8 

.98172 

1.15231 

.98172 

.95622 

.98759 

.97371 

1.08877 

.97921 

9 

-.49587 

-.60164 

-.49587 

-.47889 

-.49739 

. 00000 

.00000 

-.00000 

10 

-.96449 

-1.01358 

-.96449 

-.96440 

-.96847 

-.9737: 

-1.08877 

-.97921 

11 

-.90232 

-.89318 

-.90232 

-.90235 

-.90254 

-.92282 

-.90418 

-.92240 

12 

-.80042 

-.80067 

-.80042 

-.80041 

-.80192 

-.82844 

-.83000 

-.83013 

13 

-.66948 

-.66861 

-.66948 

-.66948 

-.67048 

-.70492 

-.70373 

-.70597 

14 

-.51252 

-.51223 

-.51252 

-.51252 

-.51328 

-.55379 

-.55348 

-.55468 

15 

-.33603 

-.33589 

-.33603 

-.33603 

-.33635 

-.38160 

-.38141 

-.38207 

16 

-.14650 

-.14649 

-.14650 

-.14650 

-.14649 

-.19468 

-.19465 

-.19478 


Table 4 (continuation) 
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C/Ax 

N 

e 

E 

e 

G 

R 

EO 

G2 

R2 

E02 

R 

R2 

0 

170 

30 

169 

77 

99 

76 

5.7 X 10"^ 

l.A X 10"^ 


Table 5. The quantities and (e = 1 x 10 for the 

initial-value problem (44.1), (49), (44.3), solved with 
both first- and second-order upwind schemes. Scheme R 
converges to the wrong weak solution. Mesh: Ax = 1/16, 
At = -^Ax. 
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■ 

R 

R2 

exact 

1 

1.00000 

.09815 

.09786 

2 

1.03596 

.28976 

.28982 

3 

1.09933 

.47012 

.47064 

4 

1.17699 

.63242 

.63337 

5 

1.25564 

.77076 

.77177 

b 

1.32417 

.87880 

.88051 


1.37431 

.95648 

.95540 


1.40069 

.98762 

.99359 


-1.4006.9 

-.98762 

-.99359 


-1.37431 

-.95648 

-.95540 


-1.32417 

-.87880 

-.88051 


-1.25564 

-.77076 

-.77177 

13 

-1.17699 

-.63242 

-.63337 

14 

-1.09933 

-.47012 

-.47064 

15 

-1.03596 

-.28976 

-.28982 

16 

-1.^000 

-.09815 

-.09786 


Table 6. Converged numerical solutions and zone- 
averaged asymptotic solution for the 
experiments of Table 5 with schemes 
R and R2. Note the expansion shock 
in the results of R. Parameters: 

K = 0, e = 1 X 10“^, At = iAx. 
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Finally, Tables 7 and 8 show the same quantitites as Tables 1 and 2 for 
experiements in which, in the first-order schemes, the source-term dependence 
of zone-boundary values was dropped. It took the schemes of Engquist-Osher 
and Godunov 15 - 25% more steps than previously to converge to a much less 
accurate solution. Roe’s scheme is now unstable; it yields an expansion 
shock where a smooth transonic transition should be; the amplitude of this 
shock grows linearly with time. The reason is that, across the zone that 
should contain the sonic point (say, zone j), we always have (see earlier 
Figure 8b) 

(50) 

SO that the scheme reduces to 

(51) Uj = u^ + nAtSj. 

8. Recommendations 

For a scalar conservation law like Burgers’ equation there seems to be 
little reason to abandon Godunov’s first-order scheme in favor of the 
first-order Engquist-Osher scheme. The slight simplification in the flux 
calculation is accompanied by a degradation of steady-shock representation 
and no significant acceleration of the convergence to a steady state. The 
simplification achieved in Roe’s first-order scheme is too drastic: the scheme 

can not be used "as is" near a sonic point. 

Both Godunov and Engquist-Osher schemes become appreciably more elaborate 
when applied to a nonlinear hyperbolic system like the one-dimensional Euler 
equations. The interpretation of the Engquist-Osher scheme as a Godunov- type 
scheme in which any shock in the solution of a local Riemann problem is replaced 
by an overturned centered compression wave, remains valid. This modification 
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- 

0 
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- 
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6.1 X 10"^ 
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CM 
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- 

.24 
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103 
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A. 7 X 10“^ 

4.7 X 10“^ 

4.7 X 10"^ 
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2 

2 

1 
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Table 7. The quantities N , E and (e = 1 x 10 for the initial- 

Co D 

value problem (44), solved with the first-order upwind schemes 
under the assumption of piecewise uniform initial values (3.1). 
(The label ”u" in Gu, Ru, EOu stands for "uniform'*.) 

Mesh: Ax = 1/16, At = JAx. 
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asymptotic solutions for the experiment 
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makes it possible, to compute ) explicitly, while Y (u ,u ) 

tiU i. K (j^ L K 

must be determined iteratively. 

The Engquist-Osher scheme, on the other hand, requires two interior 
points in a discrete stationary shock [5], while Godunov’s scheme probably 
requires only one (this has not yet been proven) . 

For the Euler equations, however, Roe’s more drastic simplification 
of Godunov’s method will pay off, even though the scheme must be somewhat 
modified in order to reject (almost) stationary expansion shocks. Spreading 
the expansion waves in the approximate Riemann solution underlying the scheme 
will have the desired effect. 

The first-order schemes, when formulated as in Section 6, have the 
potential of achieving any desired order of spatial accuracy in a steady 
numerical solution. To what degree this potential can be realized for the 
Euler equations is at present not clear. Meanwhile, the second-order two- 
step schemes, intended primarily for solving transient problems, outper- 
form the first-order schemes in obtaining a steady solution for Burgers’ 
equation. Some experiments for the Euler equations with both kinds of 
schemes are reported in [13]. 

All schemes discussed in this paper are explicit. V/hen using their numerical 

flux- functions in an implicit configuration, which may be desirable in approaching 

a steady state, the above recommendations do not automatically carry over. As 

noted by Engqulst and Osher [A], the non-smooth dependence of F^(u ,u ) on 

G L K 

its arguments in the case (iii) of a transonic shock, that is F (u ,u ) £ 

G L K 

complicates the inversion of the implicit difference equations. In contrast, 
FeqC^L’^r^^ in all cases (i) , (ii) and (iii). 
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